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Abstract 

This paper is devoted to studying Hamiltonian oscillators in 1:1:1:1 reso¬ 
nance with symmetries, which include several models of perturbed Keplerian 
systems. Normal forms are computed in Poisson and symplectic formalisms, 
by mean of invariants and Lie-transforms respectively. The first procedure 
relies on the quadratic invariants associated to the symmetries, and is car¬ 
ried out using Grobner bases. In the symplectic approach, hinging on the 
maximally superintegrable character of the isotropic oscillator, the normal 
form is computed a la Delaunay, using a generalization of those variables 
for 4-DOF systems. Due to the symmetries of the system, isolated as well 
as circles of stationary points and invariant tori should be expected. These 
solutions manifest themselves rather differently in both approaches, due to 
the constraints among the invariants versus the singularities associated to 
the Delaunay chart. 

Taking the generalized van der Waals family as a benchmark, the explicit 
expression of the Delaunay normalized Hamiltonian up to the second order 
is presented, showing that it may be extended to higher orders in a straight¬ 
forward way. The search for the relative equilibria is used for comparison 
of their main features of both treatments. The pros and cons are given in 
detail for some values of the parameter and the integrals. 

Keywords. Hamiltonian systems, isotropic oscillator, normal form, singular 
reduction, relative equilibria. 

1 Introduction 

The use of computer algebra systems for normal forms computations is considered 
at present a routine operation. As a general reference see e.g. Sanders et al. 
p9 j and Meyer et al. [26]. Nevertheless when we deal with special classes of 
differential equations, like Poisson or Hamiltonians systems which is our case, it 
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is advisable to employ specific transformations as well as tailored variables for 
those problems [26], mostly connected with the symmetries that those systems 
might possess. More precisely we are interested in perturbed isotropic oscillators 
in four dimensions (or perturbed harmonic oscillators in 1:1:1:1 resonance). For 
those Hamiltonian systems, Sanders et a 1. explain that the 1:...:1 resonance is 
one of the more complicated, due to large number of terms in the normal form. 
This proves to be a key issue in computations of higher order approximations, 
which is needed in bifurcation analysis for some values of the parameters. For this 
reason, any strategy developed to reduce the algebra involved in the normal form 
process, as well as in the subsequent analysis built on it, like relative equilibria 
and their bifurcations, is something really desirable. As we will see, we have to 
confront with systems of polynomial equations with parameters, which represent 
a real challenge, even with computer algebra assistance. 

Continuing previous work 0 E21 (231 [22 ] in the 1:1:1 resonance, and also 
mmm in relation to 1:1:1:1 resonance, we consider in R 4 xf , the symplectic 
form w = dQ A d(j and a Hamiltonian function 

'H = H 2 + EV(q,Q-,P) (1) 

where 

= 7 ^(Qi + Q\ + Q% + Q \) + 2 ^ (?i + 92 +^ 3 + qi) (2) 

defines the isotropic oscillator, with w a positive constant and £ is an small 
parameter. In the first part of the paper we simplify expressions assuming uj = 1. 
The function V is called the perturbation, where /3 is a parameter vector, which 
may include also e. Moreover we consider that systems defined by Hamiltonian 
function Eq. m have two first integrals in involution given by 

2 = q\Q 2 — Qiq 2 + Q3Q4 — Q394, 

Li = q^Qi — Q394 — qiQ 2 + Qiq 2 , ( 3 ) 

associated to which we have rotational symmetries. We use the same notation as 
in [H]. 

Although some of the methods and techniques considered in the paper may 
be applied to a large family of Hamiltonian systems defined by ( 0 ), to give details 
of those processes we focus on the uniparametric family of Hamiltonian systems 
defined by 

■Hp(Q,q)=U 2 + £H 6 , (4) 

where 

Ks(Q, q\ P) = ( q\ + q 2 +qj + qf) 

X ( P 2 {Qi +<&-(&- <&) 2 +4 (q\ + q 2 ) {q% + qV}) (5) 

where P is now a real parameter. For P = 1 we have a central potential, i.e. an 
integrable system. When the system restricts to the manifold H = 0 then it is 


2 



equivalent to the model for the hydrogen atom subject to a generalized van der 
Waals potential. For (3 = 0 this system reduces to the model for the quadratic 
Zeeman effect. When (3 = \f 1 i we have the Van der Waals system. For this reason 
we have named the system as the generalized Van der Waals 4-D oscillator. A 
search (see [j3]) for some special solutions of the Hamiltonian system defined by 
the function (J4|) reveals that there are invariant 2-tori associated to the rotational 
integrals which include straight-line orbits over the configuration space. 

Normalization and reduction using invariants and Lie-transforms are approached 
in both Poisson and symplectic formalisms. The first procedure (see m, 0 ) re¬ 
lies on the quadratic invariants associated to the symmetries, and is based on 
previous work of the authors. Enlarging previous studies on the 1:1:1 resonance 
[6,122], the Hamiltonian ([4]) is put into normal form with respect to 13.2- Consid¬ 
ering the truncated normal form a system is obtained that is invariant under the 
the § 1 -actions corresponding to 13.2 , H, and L\. These three actions together gen¬ 
erate a T 3 -action. Reduction with respect to it leads to a one-degree-of-freedom 
system. 

Hinging on the maximally superintegrable character of the isotropic oscilla¬ 
tor, the symplectic reduction is carried out a la Delaunay using a generalization 
of those variables to 4-DOF recently proposed in [21j. Based on the normalized 
equations, our studies on this system focus on relative equilibria and their evo¬ 
lution with the parameter (3. Although first order normalization (averaging) is 
quite generic, the fact that the family includes some integrable cases, makes it 
necessary to perform higher order normalization for the stability analysis of those 
cases. 

The paper is organized as follows. Section [2] summarizes part of the work 
done in m, the three steps of the Poisson reduction are presented with the 
invariant functions involved, as well as the corresponding reduced Hamiltonian 
of our model at each step. In Section [3] we proceed likewise carrying out the 
reduction, although in different order, by three symplectic transformations: Pro¬ 
jective Euler and Andoyer as well as 4-D Delaunay. They give again, but now in 
the symplectic frame, the 1-DOF thrice reduced Hamiltonian system in the open 
domain where those charts are defined. Once carried out the toral reduction by 
the symplectic variables in Section [3l Section 0] focuses on the computation of 
the normal form by Lie transform which, taking advantage of the periodic char¬ 
acter of the unperturbed flow, which allows it to solve the homological equation 
just by quadratures. Finally Section [5] gathers the analysis on relative equilibria 
which enlarges the results presented in previous papers, in particular the search 
for periodic orbits. 

2 Poisson approach. Normalization and reduction 

There is a large literature which develops both theoretical and computational 
aspects of normal forms, see for example [21 ESI [29j [26]. In this section we briefly 
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describe several normal forms for first order of system ([ 3 ]) for u = 1 thoroughly 
presented in m- In such a work, we used Grobner Bases as main computational 
tool and, as software, the computer algebra system Maple. 

A constructive geometric reduction in stages, including regular and singular 
reductions, is performed in m and only the final steps are presented here. In the 
Poisson approach the first reduction is done with respect to the H 2 symmetry. 
This is a regular reduction and the reduced phase space is homeomorphic to CP 3 
(see 0 ). Then we carry out a second reduction whose resulting orbit space is 
stratified: four dimensional isomorphic to § 2 x § 2 , and two singular strata isomor¬ 
phic to § 2 . Finally we make a third reduction, with the integral Li, reducing the 
system to a 1-DOF system on the thrice reduced phase space. Depending on the 
relative value of the integrals, they are isomorphic or homomorphic to a 2 -sphere, 
containing one or two singular points. Besides there are singular reduced phase 
spaces consisting of a single point (see Figure 2 and USD- 

2.1 Normalization with respect to the oscillator symmetry X H2 

Given the action associated to the uniparametric group defined by Xh 2 , by using 
canonical complex variables, (see m for details) it follows that the algebra of 
polynomial invariants under that action is generated by 

T; = Q 2 + (if, i = 1 , 2 , 3,4 

vr 5 = Qi Q 2 + qi 92, vr 6 = Q 1 Q 3 + 9193, TT7 = QlQi+qiqi, 

7T 8 = Q 2 Q 3 + 92 93, vr 9 = Q 2 Qi + q 2 Qi, TTio = Q 3 Qi + 93 <74, (6) 

VTll = qi Q 2 - Ql 92, 7T12 = 9i Q 3 - Ql 93, 7T13 = 9l Qi - Qi Qi, 

7+4 = Q 2 Q 3 ~ Q 2 93) VT15 = 92 <54 — Q 2 Qi, TL6 = 93 Qi ~ Q 3 Qi- 

The Xh 2 normal form up to first order in e is expressed in those invariants as 

n = n 2 + eUs ( 7 ) 

where H 2 = (tti + 7 r 2 + ^3 + 774 ) /2 = n and 

^ 6=2 [(1 — 4/3 2 )n (7 t 2 5 + 7t 2 4 + 7t 2 3 + 7t 2 2 ) 

+2 (^ 2 - l)(7rf 1 (7T 4 + vr 3 ) - 7r 2 6 (7T 3 + vr 4 )) + /3 2 n (5 H 2 - 3^) ( 8 ) 

+ 5(1 — j3 2 )n ( 7 Tg + 7 Tg + ir 2 + 7 Tg) + (/3 2 — 4)n 7 r 2 6 ] 

The reduction is now performed using the orbit map 

p n : M 8 -+ M 16 ; (q, Q) -+ (tti,--- , 7Ti 6 ) • 

The image of this map is the orbit space for the ^-action, the images of the 
level surfaces H 2 (q,Q ) = n under p n are the reduced phase spaces which are 
isomorphic to CP 3 . The normalized Hamiltonian can be expressed in invariants 
and therefore naturally lifts to a function on R 16 , which, on the reduced phase 
spaces, restricts to the reduced Hamiltonian. 
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However, in the following we will not use the invariants nf, instead, we rely 
on ( Ki,Lj , Jfc) invariants introduced in Egea [17] by the following change of co¬ 
ordinates, 


H -2 = 5(71-1 + 7T 2 + 7T 3 + 7T 4 ), 
s = 7 Ti 6 + 7Tn, 

I<1 = ^(-7Tl - 7T 2 + 7T 3 + 7T 4 ), 
K 2 = 7T S — 7T 7 , 

ET 3 = -7T 6 — 7r 9 , 

Ll = 7I"i6 — 7Tn, 

-^2 = 7Tl 2 + 7 Ti 5 , 

L 3 = 7I"i 4 — 7Ti 3 , 


j\ = \{ni - 7T 2 - 7T 3 + 7r 4 ), 

J 2 = \{ni - 7T 2 + 7T 3 - 7T 4 ), 

«/ 3 = 7T8 + 7T7, 

J 4 = 7T5 -)- 7Tio, 

J 5 = ^5 — 7T10, 

^6 — TTg 7T9) 

J 7 = 7Ti 2 - 7Ti 5 , 

J 8 = 7Ti 4 + 7Ti 3 . 


The first integrals (see Eq. O are now among the invariants defining the image. 
The first order normal form in these invariants is 
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n (5 K 2 2 + 5 K 3 2 + 2 Li 2 + L 2 2 + L 3 2 + /3 2 (5 K x 2 + L 2 2 + L 3 2 )) 
- ((4 + (3 2 ) {K 2 L 2 + K 3 L 3 ) + (2 + 3 P 2 )Kih) £ 


(9) 


The reduction of the H 2 action may now be performed through the orbit map 
PK,L,J ■■ M 8 ->■ R 16 ; (q, Q) ->■ (H 2 , ■ ■ ■ , J 8 ) . 


Note that on the orbit space we have the reduced symmetries due to the reduced 
actions given by the reduced flows of and X^. 


2.2 Toroidal reduction over CP 3 with respect to the rotational 
symmetry 5 

Let p be the S^-action generated by the Poisson flow of S over CP 3 . The functions 

H 2 , 5, Li, L 2 , L 3j K 1 , Ji 2 , iL 3 , 

are p-invariants over CP 3 . This, in turn, leads us to the orbit mapping 

p 2 : R 16 —> R 8 ; ( 711 , • • • , ttiq) —> (iLi, K 2l K 3l L 2 , L 3 , it 2 , S). 

The orbit space p 2 (CP 3 ) is defined as a six dimensional algebraic variety in R 8 
by the following relations 

K\ + K\ + K 2 + L\ + L 2 2 + Ll = H 2 2 + E 2 , K 1 L 1 + K 2 L 2 + K 3 L 3 = H 2 E (10) 

The reduced phase spaces are obtained by setting H 2 = n, 5 = £ and then the 
second reduced space is isomorphic to <S 2 + g x .S 2 ^ (see Fig. []]). When £ = 0, it 
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X 


S\ n X pt (£ = n) 



S n+i X S n-{ (0 < € < n) 
sl X Si (€ = 0 ) 

C* x S 2 n _t (-n < € < 0) 



pt X S|„ (£ = —n) 


Figure 1: Double reduced spaces S 2 + £ x A 2 ^ for different values of the integral 


corresponds to the first reduced space of Keplerian systems by the energy. The 
second reduced Hamiltonian up to first order, modulo a constant takes the form 

Ah = i [n (5 K 2 2 + 5 AV + 2 Li 2 + L 2 2 + L 3 2 + /3 2 (5 AV + L 2 2 + L 3 2 )) 

- ((4 + /? 2 ) (A 2 L 2 + A 3 L 3 ) + (2 + 3(3 2 )K 1 L 1 ) £] 

2.3 Reduction by A x = l. Thrice reduced space V n ^i 

To further reduce from S 2 + ^ x S 2 _^ to V n ^i, one divides out the 5 1 -action, p 2 , 
generated by the Poisson flow defined by L\ = ttiq — 7Tn over A 2 ^ x S 2 _^. The 
8 invariants for the L\ action on M 8 are 

M = 1 (K 2 + K 2 + L\ + A 2 ) , N = \ ( K 2 + K 2 - L 2 - L 2 ) , 

Z = A 2 A 2 + A 3 L 3 , 5 = K 2 L 3 - K 3 L 2 , (11) 

K = K\, A 2 , H, Li. 

There are 2 + 3 relations defining the third reduced phase space 

AT 2 + A? + 2M = n 2 + f, M 2 -N 2 = Z 2 + S 2 , 

KL\ + Z = n£ , V .2 = n, E = £, L\ = l . (12) 

The Poisson structure of the variables S, K and N is given in Table [I] 

For more details see m- Then, we have defined the orbit map 

p 2 :M 6 ^M 6 ;(Ad,/F 2 ,i^ 3 ,A 1 ,A 2 ,A 3 )^ (M, N, Z, S, K, L^. 
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u 

M 

N 

Z 

5 

K 

L\ 

M 

0 

AKS 

0 

—4ATV 

0 

0 

N 

-4KS 

0 

-4LiS 

—A(KM — L\Z) 

45 

0 

Z 

0 

4Li5 

0 

-ALiN 

0 

0 

S 

AKN 

A(KM - L X Z) 

4LiAT 

0 

-41V 

0 

K 

0 

-AS 

0 

41V 

0 

0 

Li 

0 

0 

0 

0 

0 

0 


Table 1: Poisson structure in (M, N, Z, S, K, L\) invariants 



Figure 2: Thrice reduced space over the space of integrals. The vector (K , N, S ) repre¬ 
sents the coordinates. The axis of symmetry of the reduced space is the K direction. 

When we fix a value of Li = l, relations m define the thrice reduced space 

V nil = { (K, S , N ) | 41V 2 + 45 2 = f(K), 

f(K ) = ((n + 0 2 ~{K + 0 2 )((n - £) 2 ~ (K - l) 2 ) } 

which is a surface of revolution, obtained by rotating \J f(K) around the axis K, 
. Thus the shape of the reduced phase space, given in Fig. [21 is determined by 
the positive part of f(K). Since 

f(K) = (K + n + £ + l)(K - n - £ + 1){K -n + £- 1){K + n - £ - Z), 
the roots are 

k\ = —l — n — £, &2 = l + n — £, &3 = l — n + £, /C 4 = — l + n + £ . 
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So f(K ) is positive (or zero) in the subsequent intervals of K: 


/ < £, -Z > £ 


k\ < k 3 < k 2 < 
k\ < k 3 < Aq < 
k 3 < k\ < k2 < 
k 3 < k\ < fej < 


fe 4 Ke[k 3 ,k 2 \ 
k 2 Ke[k 3 ,k^\ 
ki Ke[ki,k 2 ] 
k 2 Ke[ki,k 4 \ 


( 13 ) 


The Hamiltonian on the third reduced phase space is 

We.l, = ih (3/3 2 - 2) A’ 2 + 0(1 - /3 2 ) A 

+ ^(4-/S 2 )JV + » 3 (| + f) - (i 2 + f) (f + Df 

Note that the reduced phase spaces as well as the Hamiltonian are invariant 
under the discrete symmetry S —>• —S. To see how to exploit it to obtain a full 
reducction see m, m and [2¥j. We choose not to further reduce our reduced 
phase space with respect to these discrete symmetries, because the three dimen¬ 
sional picture makes it easy to access information about the reduced orbits, and 
in this way one does not introduce additional critical points (fixed points) which 
need special attention. 

In ( K , N, 5)-space the energy surfaces are parabolic cylinders. Notice that 
for /3 2 = 2/3, the function 7i is linear in the variable space ( K,N,S ). Likewise 
for (3 2 = 1 , T~L modulo constants is independent of £ and l. Moreover for f3 2 = 4, 
H is only a function of K. Since (V n £i, {•, -}3, %s,Li) is a Lie-Poisson system, the 
corresponding dynamics is given by 


^ = 2n(/3 2 -4)S, 

^ = 2[3n(3/3 2 - 2 )K + 2#(1 - /3 2 )] S, (14) 


^ = n((3 2 - 4) {K 2 - (£ 2 + l 2 + n 2 ))K 

— (3/3 2 - 2)[6nKN + 4£Z(/3 2 - 1)IV + 2Zn 2 £], 

This system can be integrated by means of elliptic functions, but after a classi¬ 
fication of the different types of flows is made, as functions of the integrals and 
the parameter. Only then we will be ready for the integration of a specific initial 
value problem. Note that in the search of relative equilibria at the thrice reduced 
level, the geometry involved is really helpful; the intersection of the Hamiltonian 
surfaces with the reduced phase space gives the trajectories of the reduced sys¬ 
tem. Then, tangency of the Hamiltonian surface with the reduced phase spaces 
provides relative equilibria that generically correspond to three dimensional tori 
in the original phase space. Details on the analysis of this system, focusing on 
relative equilibria related to (3 2 = 0, £ = l and the case of physical interest £ = 0, 
are contained in Diaz et ad. [45]. The reader should take into account that /3 2 is 
renamed in m as A. 


3 Symplectic charts for 4-D isotropic oscillators 

For our study of perturbed oscillators, we switch now from Poisson to symplectic 
formulation. More precisely, our goal is to obtain the normal form of isotropic 
oscillators by Lie transforms, in the presence of symmetries. 

Different charts are used for the symplectic treatment of the harmonics os¬ 
cillators. Among them we find complex notation [5], action-angle variables of 
Poincare type [2Sj, [23], or variations of them in the case of resonances |13j . In 
our case, as we also need to reduce by the actions associated to the rotational 
symmetries, we have made use of a symplectic chart which incorporates the func¬ 
tions ([3]) as momenta. Moreover the connection between the isotropic and Kepler 
problems, two maximally superintegrable systems, allows to introduce the Delau¬ 
nay transformation m- The result is the 4-D Delaunay transformation proposed 
by one of the authors [21) . that we will use in what follows. There is an alternative 
procedure considering the Lissajous transformation [12], which will be presented 
elsewhere. 

3.1 The family of Euler projective transformations 

Let F(p ) be a smooth real function which is positive in its domain. We consider 
the family of transformations: V£f '■ (p, <j>, —>• 92 , 93 , 94 ), dubbed as 

Projective Euler variables, given by 

41 = F (p) sin | cos —, q-i = F(p) cos | sin (15) 

<72 = F(p) sin - sill , g 4 = F{p) cos - cos , 

with (p, </>,#, i/>) € R + X [0, 2ir) X (0,7r) X (—§,§)■ For F(p) = 1, this transfor¬ 
mation defines Euler parameters as functions of Euler angles. In this paper, we 
only consider the case F(p) = ^fp. 

The canonical extension associated to the transformation (1151) is readily ob¬ 
tained as a Mathieu transformation, which satisfies ^ Qi dqi = P dp + 0 d6 + 
T dijj + <h d(f>. The relations among the momenta are given by 

P = r,(qiQi + 92Q2 + qsQs + 94Q4), 

n = (qiQi + q2Q2){ql + qp - (qzQz + q±Qp{ql + qp 

^(qi + qiM + qf) 

'F = + qiQz + q^Qi - qzQP, (16) 

^ = 2 (qzQi ~ qiQ2 + q±Qz — qzQp, 
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Later on we will need the inverse transformation given by 


2 , 2 i 2 , 2 

P = Ql + 9 2 + 9 3 + 9 4 > 


sin# = 
sini/> = 


sin <f> = 


2 V(q( + qlM + qt) 

9? + 9 2 2 + $ + 94 ’ 

9 i 93 + g 2 g 4 

'/]'/:i - 9294 

V / 5f+«IMTS)’ 


COS 0 = 

cos if = 


COS (f = 


93 + 94 - 91 - 9 2 

9 ? + ql + Ql + qV 

9194 ~ 9293 

VW + ^Kii + a’ 

9194 + 9i4i 
V(9?+9 2 2 )(9l + 9l)' 


(17) 


Proposition 3.1 T/ie isotropic oscillator "H 2 is expressed by mean of the Euler 
projective variables as 


n^p, 6, -, -,P, *, 0 , $) = ^ + 2pP 2 + 


0 2 + 


vj/2 _|_ $2 _ 2 $ \p cos 0 


r.2 a 


(18) 

Proof : Considering any given initial condition for the Hamiltonian Hu = h, by 
applying relations (fT5l) and (HI, we obtain, by straight forward computations, 
the above expression for the isotropic oscillator % 2 . q.e.d. 


Note that the Hamiltonian has two cyclic variables, which manifests the two 
symmetries associated with our system we have refer in Q 


3.2 From Euler to Projective Andoyer variables 

Andoyer mi symplectic variables are well known in rotational dynamics and re¬ 
cently in attitude and control. Most often they are denoted by (/3 2 , p, v, A, M, N ) 
or (h, g,£, H,G, L). In the following, in order to avoid confusion with invariants 
notation of previous section, we propose to use (iq, 112 , 1 / 3 , Hi,H 2 ,H 3 ) for refer¬ 
ring to them. Moreover, as in 3-D, it is convenient to introduce the following 
functions 

Ci = coscji = U 1 /U 2 , c 2 = coscr 2 = H 3 /H 2 , (19) 

with Si = \J 1 — c 2 . Hence the transformation from Euler to Andoyer (<^>, 9, if, 
$,e,^) -> {ui,U 2 ,u 3 ,Ui,U 2 ,U 3 ) is given by 


COS 9 = ClC 2 + Sl s 2 cos iz 2 , 


sin(^ — iti) 


sin U 2 
sin# 


« 2 , 


sin (<f - u 3 ) 


sin U 2 
sin# 


si, 


<h = C/ 3 , 
* = U!, 


( 20 ) 


Q = U 2 \ 1 - 


c 2 + c| — 2 ci c 2 cos ( 


sin 2 6 


with |Hi| < U 2 and |H 3 | < H 2 . By adding the variables (p,P) we get a 4- 
D set (p, ui, U 2 , u 3 , P, Hi, H 2 , H 3 ) of symplectic variables that we call ‘Projective 
Andoyer variables’. 


10 



















Theorem 3.2 In Projective Andoyer variables, the system defined by (0), reg¬ 
ularized by ds = 1/(4 p)dr, includes the Keplerian system for any value of the 
integral U 3 . 

Proof : Following Poincare technique, we first regularize m and our new 
Hamiltonian will be K, = (PL^ — h)/(4p). Considering now Proposition 13.11 after 
some straightforward calculations, the Hamiltonian of the 4-D isotropic oscillator 
is given by 

(2i » 

in the manifold A = —oj/8 and 7 = j. But the above Hamiltonian corresponds 
to the Kepler system in polar nodal variables, see |10j . q.e.d. 

3.3 Delaunay symplectic chart 

Following m we plan to normalize perturbed isotropic systems by Lie-transforms 
a la Delaunay in the following Section. Thus, we still need to implement (in part 
of the phase space) another canonical transformation of Hamilton-Jacobi type: 
D 7 : ( L,G,£,g ) —>• (p, U 2 , P, H 2 ). The process hinges on the H-J equation built 
on the function (1211) . For details of this transformation T > 7 we refer to m here 
we only include the final expressions. Among the relations defining Delaunay 
transformation, which can only be given in implicit form, we take 


U2 = 9 + /, 

II 

£ 


p = a (1 — e cos E), 

Le sin E 
a (1 — e sinK) ’ 

( 22 ) 

where a = a(L), e = e(L, G ) are given by 



« = l 2 / 7 , V = G/L, 

e = y 7 ! - 

(23) 


and / and E are auxiliary angles. The symplectic variable l is related to them 
by 

£ = E — esinE, tan//2 = y/(l + e)/(l — e) tanK/2. (24) 

Later on we also need 1/p = (1 — e cos f)/p, where p = arj 2 , and 

p cos / = a(cos E — e), p sin / = a/q sin E, d£ = (1 — e cos E) dE. (25) 

Completing the functions of the momenta given above (1231) . it is convenient to 
introduce two more state functions 

w = Ui/L , z = U 3 /L. (26) 

The previous results may be summarized in the following theorem: 
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Theorem 3.3 By composing the three symplectic transformations 


/ qi, 

<? 2 , 

Q3, 

QA \ Projective Euler 

( 

<t>, 

e, 

\ 

\Qi, 

Q2, 

Q3, 

Q4J 

v P, 


0 , 

'F ) 


PA 


( 

Ul, 

g, 

U 3 \ 

Delaunay 

( p 

Ul, 

U'2, 

U 3 \ 

\ L, 

Ui, 

G, 

u 3 ) 


\ p, 

Ui, 

U 2 , 

U 3 ) 


together with the regularization ds = l/(4p) dr, the Hamiltonian of the 4-D 
isotropic oscillator is given by 


Ho — 



(27) 


in the manifold Ho = — ui/8. 


The set of variables (£, g, ui,u 3 , L, G , XJ\ , C/ 3 ) is what we call 4-D Delaunay chart. 
As in 3-D with the classical Delaunay variables, they represent a generalized set 
of action-angle variables (see [28]) • In the next section we show the interest of 
this symplectic transformation. For the benefit of the reader let us mention that 
Moser and Zehnder m introduced action-angle variables for the Kepler problem 
in M n and called them Delaunay variables. When we restrict to n = 4, those 
variables do not coincide with the set built in this paper. 

Finally, in order to compare results of our analysis of relative equilibria, we 
need to establish the connection between the variables K (Poisson approach) and 
G (symplectic approach) defining the thrice reduced space. After some algebraic 
manipulations, we obtain 


4 G 2 = ^(n 2 + e 2 + / 2 )- \k 2 ~N, 


(28) 


which expresses G as function of K and N, invariants defining the thrice reduced 
space, and the first integrals. 


4 Delaunay normalization of perturbed 4-D isotropic 
oscillators 

Our goal in this section is to make normalization by Lie transforms [9], using 
4-D Delaunay variables. As we will see they allow to express perturbed isotropic 
oscillators in a ‘convenient’ form, to implement Lie transforms up to higher order 
in an efficient way. What we mean by convenient is the following. 

It is well known that the isotropic oscillator Ho, as the Kepler system, defines 
a maximally superintegrable Hamiltonian system (see Fasso [20]), whose flow 
is made of periodic orbits. In that case Cushman [3] proved that any smooth 
function F over its phase space may be decomposed in a unique way into a sum 
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F = F° + F* with the following properties: (i) { F °, "Ho} = 0. In other words, F° 
belongs to the kernel of the Lie derivative generated by Ho'- £o : F —> {F, 'Ho}', 
(ii) There exists a smooth function F such that {T 1 , HLo] = F*; in other words 
F* belongs to the image of the operator Cq. In particular when the Lie transform 
(see [9], [26]) 

£ : [l^uu^L&Uulh) -7 U[,U’ 3 ) 


is carried out a la Delaunay m, the homological equation, which relates the new 
Hamiltonian /j\)T-L'j with the generating function ^{e 3 /j\)W}, is given by 

{W ] ,'Uo} + 'H' j =U j , j> 1. (29) 

Then, considering the splitting coming from Cushman theorem, the equation (1291) 
may be solved choosing 

1 f 27T ~ \ f ~ 

n'j = <Uj>= — J Uj d£, Wj = -J (Uj- <Uj>) dl, ( 30 ) 

where, according with (|27|). dH.o/dL = y 2 /L 3 . The Hamiltonian ([TD of perturbed 
isotropic oscillators in Delaunay variables is given by 

7 2 

n{l,g,ui,u 3 ,L,G,Ui,U 3 ) = + £V{l,g,ui,u 3 ,L,G,Ui,U 3 ', j3). 

and the normalized Hamiltonian up to order k takes the form 


%' = - 


7 


2 L' 2 


+ ^-Hj(-,9',u , 1 ,u' 3 ,L',G', u i u ! i] p) + 0(e k+1 ). 


The main feature of the process is that, at each order, solving the homological 
equation (129[) only involves quadratures. 


4.1 Delaunay normalization and symmetries 

Apart from the procedure associated to Delaunay normalization we have just 
mentioned, when we restrict to systems with the symmetries Q, the use of 
Delaunay chart shows its full advantage. Indeed, in that case u\ and u 3 are 
cyclic, in other words, we have 

m, 9 , L, G, U 1} U 3 ) = + eV(l, g, -, —,L, G, U \, U 3 ;fF, e). (31) 

In geometric language, the use of Delaunay variables has carried out the toral 
reduction associated to the actions defined by the symmetries ([3]) , and the func¬ 
tion (1311) is the reduced Hamiltonian. The Hamiltonian system given by (1311) is 
a 2-DOF system. Thus, the normalized system will take the form 


W = - 


7 


k 

+ Y, a', - -,L', G', U[, t/'; P) + Q(e k+1 ). (32) 


2 L' 2 


n | 
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In our case normalizing up to the order needed, after truncating, we will 
obtain an integrable 1-DOF Hamiltonian system. From now on we drop primes 
in the variables in order to simplify the expressions. 

4.2 Carrying out the normalization. The van der Waals model 
as a benchmark 

In order to implement these normalizations, since we are dealing with polynomial 
perturbations, from experience with 3-DOF Kepler systems we introduce the 
auxiliary variable E by using the following proposition. 

Proposition 4.1 The functions F(p,6) = p m cos n 6, ( m,n € N and m > n), 
expressed in Delaunay variables, could be written as F(p,9) = YaL 0 (CiC°si E + 
Si sin i E), where Ci and Si are given by expressions of the form Yj( c ij cos j g + 
Sij sin j g) which belong to the kernel of the Keplerian Lie derivative, and Cij, Sij 
are rational functions of the momenta. 

Proof : It is straightforward, based on relations given in (l20l) . (|22l) and (l25l) . 

q.e.d. 

Applying the above result, the Hamiltonian (j4j) takes the form 
2 fc j 

U = -^2 ! Y.~ Y. {C L cosiE + S H s[niE ) + 0{e k+1 ), ( 33 ) 

j =0 J ' i 

where functions Cij = Cij(g, L, G,Ui,U 3 \ (3) and Sij = Sij(g, L,G,U\,U 3 ] /3) be¬ 
long to the kernel of the Lie transform. 

Thus, by considering the differential relation (1251) and by replacing in (1301) . 
we obtain the part in the kernel as follows 

1 f 27r 

<TLi>=— / Thdl 
Jo 

= — / (Cn cos iE + Sn sin iE)(l — e cos E) dE 

2-k Li ^' 

J{) »=o 

= Coi + Cn cos g + C 21 cos 2 g. 
with Ci 1 = Cn(L, G, Ui,U 3 -/3) given by 

Coi = —« 2 ( 2 + 3e 2 )[2 + (/ 3 2 — l)(2cfc| + s 2 s 2 ) ], 

Cn = -^a 2 (/3 2 - l)(4 + e 2 )ecic 2 sis 2 , (34) 

C 21 = t^« 2 (/3 2 - 1 )e 2 44- 

lb 

Then, we solve the first order homological equation (1291) . The first order normal¬ 
ized Hamiltonian takes the form 
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T~L — — ^2 + e (Coi + Cn co sg + C 21 cos 2 g) + 0(e 2 ) 

The first order W 1 of the generating function is obtained computing the second 
quadrature in (l30l) . then, taking into account the shorthand a = (3 2 — 1, we get 
the following expression 

Wi = ( a ( 2c i c 2 + s2s i) + 2 ) 

^ci (e 3 sin(3'«) — 9e 2 sin(2 E) + 3e(3g 2 + 5) sin(it)) + as\S 2 

+siS 2 (l — v) 2 { ~ 15esin(2g — E) + (9 + 6 rj) sin(2g — 2 E) 

—esin(2g — 3 Ej) 

— siS 2 (l + r/) 2 (l 5 esin (25 + E) — (9 — 6 rj) sin(2g + 2 E) 

+esin(2 g + 3 Ej) 

+ 4 ciC 2(1 — 77 ) (e 2 sin (g — 3 E) — 3e(3 + rj) sin (g — 2 E) 

+3(2 rf + rj 2 + 5) sin (g - E )) 

— 4 ciC 2(1 + 77 ) (e 2 sin (g + 3 E) — 3e(3 — rj) sin(g + 2 E) 

+3(-2 if + g 2 + 5) sin (5 + E))^ . 

Analogously we obtain the second order normalization for (1321) . which will be 
used in the study of stability. The second order coefficients are 

/'“y _ 2 2 2 2 2 6 

G 02 = —60a C 2 c 1 s 1 s 2 e 
2 

+ (^-(1139c 2 4 c? — 1486(c 2 4 c 2 + c 2 2 cf) + 579(c 4 + c 2 4 ) 

+1268 C 2 2 c\ - 518(c? + 2 c 2 2 ) - 61) + 63 (a A 3 ,i - 1)) e 4 
+ (a 2 (2089 c 2 4 c 4 - 1756(c 2 2 c? + c 2 4 c \) + 399(c? + c 2 4 ) 

+2048 C 2 2 c\ - 628(0? + c 2 2 ) + 229) - 396 (a A 3a - 1) ) e 2 
+ (a 2 (557 c 2 4 c? - 382(c 2 4 c? + 382 c 2 2 cf) + 17 (cf + c 2 4 ) 

+308 c 2 2 c 2 - 22(c 2 2 + cf) + 5) - 96 (a A 3)1 - 1) ), 
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Ci 2 = C 2 C 1 S 1 S 2 (( - 155a 2 (A 2j i + 41cfc| + 36) + 192 )e 5 

-(447a 2 (A 2j1 - 19 c\c\ - 290) + 708a + 2)e 3 
- (67 a 2 (A 2j1 + c\c 2 2 + 27) - 1200 a - 16) e), 

C22 = sis 2 (^ - 60a 2 c 2 2 c 2 e 6 

+ (—^- (71A 5>2 + 8c 2 c| — 14) — 57a )e 4 
+ (a 2 (— 147A 2j1 + 25c 2 C 2 - 66) +486 a + 3)e 2 ), 

C 3 2 = cisfc^a 2 ^85e 5 + 10e 3 ^ , 

^ 95 2 2 2 4 

C42 = a z s{s z 2 e A 

8 

where 

A n ,i = cf + c|(—l)*?rc 2 c|(—1)7 

As the van der Waals system includes several integrable cases, in order to 
study stability in those cases, higher order normalization, at least second order, 
is needed. This subject will be tackled in future works. 

Here we try only to identify relative equilibria. In fact, this is the procedure 
to follow in order to classify and compute different types of orbits defined by 
our system. For the generic case, first order normalization is sufficient. The 
second order normalization is required when we consider degenerate situations, 
connected with integrable cases. For this reason we give both, first and second 
order in Section [4.21 

5 Searching for relative equilibria and invariant tori 
of the first order normalized system 

Focusing on the task that we have announced before, we use the first order 
normalization of the Hamiltonian. Thus the corresponding differential system is 
given by 

(35) 

(36) 

(37) 

(38) 
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dU 

7 2 

dV 

L = 

dn 

= 0 , 

£ = 

dL ~ 

h 

+ £ 5l’ 

d£ ~ 


dn 


dV 

n 

dn 


9 - 

dG 


e dG ’ 

Lt — 

dg ’ 



dU 


dV 

Ui = 

dn 

= 0 

u\ ■ 

dUx 

— 


du\ 


dU 


dV 

u 3 = 

dn 

= 0 

U 3 ■ 

' dU 3 

— 

£ du 3 5 

du 3 







In other words L, U\ and U 3 are integrals, as we already know, and the system 
splits into a 1-DOF Hamiltonian system, namely 


dV 

9 ~ £ dG 1 

and three quadratures coming from 



&P_ 

“ 13 +e aL ; 


U\ = £ 


dV 

Wi 


U 3 = £ 


dV 

dU 3 


(39) 


as soon as the system Eqs. (1391) is solved. We are not going to explore this way, 
which moreover involves hyperelliptic integrals. In other words, it will be more 
convenient to rely on numerical integrations. 

As we have already said, our previous work on this model concentrated on 
searching for relative equilibria related to singular points of the energy-moment 
map. We have already mentioned that the complete analysis should be done 
using invariants. In fact the lower dimensional relative equilibria, i.e. those that 
correspond to invariant S 1 or T 2 , are given by the singularity of the moment map 
for the T 3 -symmetry group, and can be described by a moment polytope; for 
details see [ 1 TI| . 

Excluding those solutions, in the remaining open domain we may use sym- 
plectic charts. Relative equilibria of the system (1351) - (1381) may be classified as 
follows: 

• invariant 3-tori. They are the solutions (g, G ) of the system defined by 


G 



■ dV n 
a = e 8G =0 - 


(40) 


as functions of Hi, U 3 and (5. 

• invariant 2-tori. In a similar way, the search for invariant 2-tori splits into 
the search the roots of the subsystems. Namely 


dV 

G = -£— = 0 , 

og 


• dV n 
9 = £ dG =°> 


dV n 
u l = £ 7 — = 0 , 


'dlh 


where the possible roots will be function of U 3 and j3, and 


dV 

G = -£— = 0 , 

og 


• dV n 
9 = £ dG =°> 


U 3 =£ 


dV 

dlh 


= 0, 


(41) 


(42) 


where the possible roots will be function of U\ and /3. 

• periodic orbits: The research for periodic orbits of the original system is 
equivalent m to find the roots of the first reduced system of equations 

dV dV dV dV 

G = -e-z- = °> 9 = = 0, *1 = = 0, u 3 = e-v- = 0, (43) 

og oGr oU\ 0U3 

Solutions of the system given by Eqs. (SOD define invariant 3-tori T 3 (£,ui,u 3 ), 
except for the singular points m- The system defined by Eqs. m gives the 
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invariant 2-tori T 2 (£, u^), each of them is attached to a 2-torus T 2 (g, u\) made 
of fixed points. Likewise the system defined by Eqs. (I42|) is the invariant 2-tori 
T 2 (£, u\), each of them is attached to a 2-torus T 2 (< 7 , 1 * 3 ) made of fixed points. 
Finally, the intersection of the previous types of 2-tori, defined by the system 
of Eqs. (1431) . correspond to periodic solutions S x (f), each of them is attached 
to a 3-torus T 3 (< 7 , zii, 7 / 3 ) made of fixed points. They correspond to the short 
period solutions related to the unperturbed system. In what follows we search 
for invariant 3-tori and periodic solutions of short period. 


5.1 Searching for invariant 3-tori 

To find families of relative equilibria (periodic orbits, invariant tori, etc), we start 
searching for invariant 3-tori, the common condition for all the cases. Explicitly, 
according to above paragraphs, the equations (T4U1) are 


dC 01 

dG 


dC u 

a g COS9 + 


dC 2 1 
dG 


cos 2 g = 0 , 


(44) 


(Cn + 4 C 21 cos g) sin g = 0. 


(45) 


Due to the structure of Eq. (|45l) . and keeping in mind the domain of existence 
of 4-D Delaunay variables, the only possibility leading to roots is sin g = 0. 
Replacing in (1441) . we have an equation to be solved for each of the values of 
cos g, namely ±1. Such equations are: 

8C 01 dCu dC 2 1 8C 01 dCu dC 2 1 

dG + dG + dG ’ dG dG + dG ' 1 ’ 

Both equations can be solved at one time applying the following strategy. By 
using the state functions (l23l) and (1261) . both equations (jj6|) can be rewritten 
in the form; R(g, z,w)x(r], z,w) + Q(rj,z,w) = 0 and R(g, z,w)x(rj, z,w) — 
Q(rj, z, w) = 0, in which: 

x(rj,z,w) = \/(l — w 2 )( 1 — z 2 )e, 

R(rj, z, w) = (—(3 + 4 a)g 6 + (5w 2 + 7 w 2 z 2 + 5 z 2 ) arj 2 — 20 w 2 z 2 a) rf, 

Q(r 1 , z,w ) = a w z [rf — (10 + llu> 2 + w 2 z 2 )i] 4 — 11 z 3 if (47) 

+ (15u; 2 + 17 w 2 z 2 + 15 z 2 )if + 5 z 5 rj — 20u; 2 2: 2 ] . 


By multipliying equations Rx+Q = 0 and Rx—Q = 0 we obtain a new polynomial 
P(rj) = R 2 x 2 — Q 2 in g as independient variable, where w, z and a are taken 
as parameters. The roots of P represent the set of solutions for the original 
equations (1461) . Explicitly the polynomial is: 

P{g) = a 6 ry 6 + a 5 g 5 + a 4 ? f + a 3 ? f + a 2 g 2 + a\g + a 0 , (48) 
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where the coefficients are given by 
dg = — 9 — 24 a — 16 a 2 , 

«5 = tc 2 (16a 2 + 9 — a 2 z 2 + 24a) + 9 + 24a + 16a 2 + z 2 (24a + 9 + 16a 2 ), 

04 = w 2 ( 24a 2 + 18az 2 + 6 a —9 —9 z 2 + 30 a 2 z 2 w 2 ) + z 2 ( 6 a + 24a 2 — 9), 

03 = w 4 {—40 a 2 — 42 a z 2 — 30 a — 34 a 2 z 2 + 2 z 4 ) + 9 z 2 w 2 — 30 z 2 — 40 a 2 z 2 

— 40 a 2 z 4 + arc 2 (198 z 2 — 34az 4 + 30 + 285 az 2 + 40 a + 42 z 4 ) — 30 a z 4 , 
02 = a (30 w 4 + 42 z 4 ie 4 + 192 z 2 w 4 + 15 a re 4 — 17 a z 4 rc 4 + 266 a z 2 u > 4 — 30 z 4 

+ 266 a z 4 rc 2 + 192 z 4 rc 2 + 180 z 2 w 2 + 290 a z 2 w 2 + 15 a z 4 ), 
ai = a 2 (25 w 4 - w 6 z 6 + 27 w 6 z 4 - 51 ic 6 z 2 + 25 w 6 + 27 z 6 rc 4 - 139 z 4 w 4 

- 225 z 2 w 4 + 25 z 6 + 25 z 4 - 51 zV - 225 z 4 w 2 - 50 z 2 w 2 ) 

+ 150 z 2 w 4 + 150 z 4 w 2 + 162 z 4 w 4 , 

oo = 5 a (a ic 6 z 4 — 3 a rc 6 z 6 — 5 a w 6 + 7 a ic 6 z 2 + 24 z 4 w 4 + a z e w 4 

+ 5az 2 rc 4 + 18az 4 rc 4 — 5az 6 + 5 az 4 w 2 + 7az 6 rc 2 ). 

This polynomial corresponds with the one obtained from the system (|14[i when 
we look for relative equilibria within the Poisson approach m ■ The complexity 
in both analyses is similar because the polynomials are of the same degree. It is 
straighforward to check that when we impose the constraint z = 0 , we recover 
the expressions of the study done with the classical Delaunay variables jl9j . 

The search for possible invariant T 2 solutions of Systems (14ip and (14211 requires 
a similar analysis. We will not refer to them here. Instead, we focus on periodic 
orbits, where we obtain a benefit working in symplectic formalism. 


5.2 Searching for periodic orbits 


With Poisson formalism, finding periodic solutions is in correspondence with com¬ 
puting relative equilibria of the system defined with the normalized Hamiltonian 
([9]). This study requires to deal with the constraints defining CP 3 and leads 
to a large polynomial system; for details see m- Although partial results are 
obtained, the presence of the physical parameter introduces complicated expres¬ 
sions. 

The situation changes dramatically when we switch to symplectic formalism. 
The complete system (1431) related to periodic solutions takes the form 


dC 01 

dG 

dCo 1 
dUx 
dC 0 1 
dU 3 


dC n dC 2 1 

+ ~dG~ COS 9 + ~dG~ C ° s29 = °’ 
(C 11 + 4 C 21 cos g ) sin g = 0, 
, 9 C n , dC 2 1 o n 

+ W COS5 + W COS25 = °’ 
+ ~W ~ COS9 + ~w~ cos2a = °- 

dU 3 dU 3 
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The solutions of the system are the set of the periodic orbits in the first order 
normalized problem. By solving the above system imposing a / 0, we obtain 
that periodic orbits are characterized by U\ = O = 0 and 


I jj | = I jj | = / e cos ff (4 + 5 e cos g + e 2 ) + 1 — e 2 
1 3 y ecosg (5 ecosg + 8 + 2e 2 ) + 3 + 2 e 2 

I T j | = lTJ | = / ecos g (-4 + 5 e cos g - e 2 ) + 1 - e 2 

1 3 Y ecosg (5ecosg — 8 — 2e 2 ) + 3 + 2e 2 

Then (14311 becomes a system of three equations given by 

dc m , dCu , 96 2i 0 n 

^ + ^co 59 +^c°s2 s = 0, 

(Cn + 4 C 21 cos g) sin g = 0, 

<90,1 , r9C n , dC 21 0 n 
+ cos <7 + cos = 0. 


(49) 


(50) 


90 


dU, 


90 


Again we look for solutions when sin <7 = 0. Then, we drop the second equation 
and we have two systems to solve related to cos g = ±1: 


and 


w 


(**) 


96 01 , 

9C n ( 

9C 2 i n 

9C 0 i , 

967n , 

96,, 

+ 

dG 

+ 

dG 

= 0 

9G ’ 

+ 

90 

+ 

90 

90 

96 01 

dC n , 

dC 2 1 n 

9C 0 i 

96/, ( 

96 21 


dG dG dG ’ 90 90 90 

Considering case (z), from the second equation we obtain 

c 2 (e 2 + 3 e + 1 — c| (2 e + 3) (e + 1)) =0. 


= 0, 


= 0. 


(51) 


(52) 


Thus, we have either c 2 = 0 or 

2 1 + 3 e + e 2 

C2 “ (3 + 2 e) (1 + e)' 


(53) 


If c 2 = 0, we obtain a = —3/4. Otherwise, from the first equation we find 


3e (3 + 2 e) 2 

3 e 4 + 2 e 3 — 10 e 2 — 3 e + 8 


(54) 


Likewise, for the case ( ii ) the equation is 

c 2 (e 2 — 3 e + 1 — c 2 (2 e — 3) (e — 1)) = 0. 
Hence, again c 2 = 0 or 


1 — 3 e + e 2 
(2 e — 3) (e — 1) ' 


(55) 
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If C 2 = 0, we obtain a = —3/4. Otherwise, from the first equation now we find 

3e (3 — 2 e) 2 

3 e 4 — 2 e 3 — 10 e 2 + 3 e + 8' (56) 

Thus, for each member of the Van der Waals family defined by a, we have the 
relation between integrals which give periodic orbits. The period is computed 
replacing those values in the right hand terms of (13511 . 


Conclusions and future work 

In contrast to some claims in [5], Poisson and symplectic formulations of per¬ 
turbed 4-D isotropic oscillators systems behave as complementary in the analysis 
of the main features of those systems. The first approach, built on the invariants 
and their algebraic relations, is necessary in order to carry out regular and sin¬ 
gular reductions of the systems, in particular when singular points are present in 
the reduced spaces. Nevertheless the use of those invariants and their relations, 
introduce large computations in the search of relative equilibria, specially when 
physical parameters are involved. Only those equilibria related to the symmetry 
groups of the system are more easily computable. 

The situation is rather different when the symplectic approach is followed. 
From the geometric mechanics point of view this helps to portrait the toral 
structure of the phase space, necessary when stability KAM theory is applied. 
Moreover, from the algebraic perspective, the number of equations and relations 
involved falls sharply because there are no constraints. Normal forms are com¬ 
puted very efficiently in this frame. Of course, this is at the expense of considering 
only an open domain of the phase space, where those symplectic variables are 
defined. For this reason we should begin with Poisson formalism. It is only af¬ 
ter studding singular points that we can switch to symplectic techniques with 
properly chosen variables. 

As an illustration, the Van der Waals family is studied in detail here, showing 
the pros and cons of both approaches. The reconstruction process connecting the 
whole analysis of relative equilibria with the original system, is still to be tackled. 
In particular, the analysis of the different types periodic orbits. 

An open question is whether there are other invariants which can lead to 
similar equations to those appearing in symplectic variables . In the same vein, 
connected with our choice as 4D-Delaunay variables [23], it might be worth con¬ 
sidering other recent proposals [17], [30]. Results on this, now in progress, will 
be published in [3]. 
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